A FAST ALGORITHM FOR MACMAHON'S PARTITION ANALYSIS 



GUOCE XIN 

Abstract. This paper deals with evaluating constant terms of a special class of rational 
functions, the Elliott-rational functions. The constant term of such a function can be read off 
immediately from its partial fraction decomposition. We combine the theory of iterated Lau- 
rent series and a new algorithm for partial fraction decompositions to obtain a fast algorithm 
for MacMahon's Omega calculus, which (partially) avoids the "run-time explosion" problem 
when eliminating several variables. We discuss the efficiency of our algorithm by investigating 
problems studied by Andrews and his coauthors; our running time is much less than that of 
their Omega package. 



1. Introduction 

Zeilberger [20] proved a conjecture of Chan et al. [11] by proving an identity equivalent to 
(1-1) CT---CT -— — - = Ci---Cn~i, 

where C^s are the Catalan numbers. 

This identity should be interpreted as taking iterated constant terms [10]; i.e., in applying 
CT^^ to the displayed rational function, we expand it as a Laurent series in the result is 
still a rational function and we can apply CTx„_i, • • • , CT^.^ to it iteratively. 

The idea behind the above treatment is to give a proper series expansion of l/(xj — xj) for 
every i and j, so that all of the expansions are compatible. Once we have determined the 
relations between the x's, there is no confusion about their series expansion. For instance, we 
can let 1 > Xi > ■ ■ • > x„. For the particular rational function in equation (1.1), which is 
symmetric in the x's, there is no confusion after a total ordering on the a;'s is given. 

Here we present a slightly different, but more efficient, approach, by means of applying the 
theory of the field of iterated Laurent series. We first treat the rational function in question 
as an iterated Laurent series, by which we mean we expand it as a Laurent series in then a 
Laurent series in and so on. Then we take the constant term. This idea led to the study 
of the field of iterated Laurent series in [18, Ch. 2], which applies to MacMahon's Partition 
Analysis. 
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MacMahon's Partition Analysis is suited for solving problems of counting solutions to lin- 
ear Diophantine equations and inequalities. Using MacMahon's approach, problems such as 
counting lattice points in a convex polytope, counting integral solutions to a system of linear 
Diophantine equations, and computing Ehrhart quasi-polynomials, become evaluations of the 
constant term of an Elliott-rational function: a rational function whose denominator has only 
factors of the form A — B, where A and B are both monomials. An example of such is the 
rational function in (1.1). 

MacMahon's technique has been restudied by Andrews et. al. using computer algebra in 
a series papers [1-9]. New algorithms have been found and computer programs such as the 
Omega package have been developed. 

The constant term (in one variable) of an Elliott-rational function can be read off imme- 
diately if its partial fraction decomposition is given. However, the coefficients of a rational 
function must lie in a field to guarantee the existence of its partial fraction decompositions, and 
the classical algorithm for partial fraction decomposition is rather slow because the coefficients 
contain many other variables. The above two problems are solved by applying the theory of 
iterated Laurent series and a new algorithm for partial fraction decompositions in [17]. 

In section 2, we give the basic theory of iterated Laurent series. The fundamental structure 
theorem tells us when a formal Laurent series is an iterated Laurent series. In section 3, we 
introduce MacMahon's partition analysis. In section 4, we develop an efficient algorithm for 
MacMahon's partition analysis by combining the theory of iterated Laurent series and a new 
algorithm for partial fraction decompositions. The theory of iterated Laurent series is crucial 
in avoiding the "run-time explosion" problem [7, p. 9] when eliminating several variables. In 
section 5, we use our Maple package to test the efficiency of our algorithm. We investigate 
problems related to fc-gons, generalized Putnam problems, and magic squares [4; 6; 7; 9]. The 
known formulas are obtained within seconds, and several new formulas are produced in minutes. 
Finally in section 6, we point out several ways to accelerate the computer program. There are 
also ways to make the computation easier that are hard to implement on the computer. As an 
example, we give a simple proof of the formula for /c-gon partitions in [4] . 

2. The Field of Iterated Laurent Series 
By a formal Laurent series in xi, . . . , x^, we mean a series that can be written in the form 

oo oo 
oo in=~oo 
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where aii...i„ are elements in a field K. For formal Laurent series, the definition of the constant 
term operator is clear: 

Definition 2.1 (Natural Definition). The operator CT^;^. acts on a formal series m. xi, . . . ,Xn 
with coefficients ai^^,,,^i^ in K by 

^ (n,...,jn)GZ" (n,...,i„)eZ",jj=o 

The simplest way to apply the natural definition would be to work with all formal series 
X^jj i aii,...,j„a^i^ ■ ■ 'X'^ni where . . . ,in) ranges over all elements of Z". Unfortunately, they 
do not form a ring. Therefore we usually work in a ring, such as the ring of Laurent series 
K{{xi, . . . , Xn))'- formal series of monomials where the exponents of the variables are bounded 
from below. But wc need a larger ring or even a field that includes all rational functions, 
because many constant term evaluation problems involves rational functions. 

Let X be a field. We define K{{xi)) to be the field of Laurent series K{{xi)), and define the 
field of iterated Laurent series K{{xi, . . . , inductively to be K{{xi, . . . , Xn-i)){{xn)), which is 
the field of Laurent series in x„ with coefficients in K{{xi, . . . , Thus an iterated Laurent 

series is first regarded as a Laurent series in Xn, then a Laurent series in Xn-i, and so on. An 
iterated Laurent series obviously has a unique formal Laurent series expansion. However, it 
is not obvious which formal series are in K{{xi, . . . ,Xn))- The fundamental structure theorem 
solves this problem nicely. 

Wc define a total ordering ^ on monomials by representing x]^ ■ ■ -x^^ by (ii, . . . G Z", 
where Z" is ordered reverse lexicographically. So x\ -< Xj for all i < j and s G Z. We define 
the support of a formal Laurent series by 

supp (^iu.-,inXi •••<":={ (ii, ... in) I aiu-,in 7^ }. 

(n,...,i„)eZ" 

Recall that a totally ordered set S is well-ordered if each nonempty subset of S contains a 
minimal element. 

Theorem 2.2 (Fundamental Structure). A formal series in xi, . . . ,Xn belongs to 
K{{xi, . . . , Xn)) if and only if it has a well-ordered support. 

The proof of this theorem is omitted. For details, see [18, Proposition 2.1.2]. The result 
gives us an overview about when a formal Laurent series is an iterated Laurent series. 
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The fundamental structure theorem, together with the simple and useful fact that any 
subset of a well-ordered set is well-ordered, justify the application of the natural definition in 
K{{xi, .... ,T„)) because of the following three properties: 

PI. GTx^ ■ K{{xi, . . . , Xn)) K{{xi, . . . ,Xi, . . . , This property is necessary to make 

the natural definition applicable. 
P2. CTj;^ — J2i ^^xk ^i- This property is the key to converting many problems into 

simple algebraic computations. 
P3. CTj;. CTxj F — CTxj CTj;. F. This property may significantly simplify the constant 

term evaluations. 

We define the order ord(/) of an iterated Laurent series / to be the minimum of its sup- 
port, which is well-ordered by the fundamental structure theorem. We have the following 
composition law. 

Proposition 2.3 (Composition Law). Suppose that f belongs to K{{xi, . . . , Xn)) andord{f) > 
ord(l). Then for any hi ^ K for all i, 



is well defined and belong to K{{xi, . . . ,Xn)), in the sense that all of its coefficients are finite 
sum of nonzero elements in K. 

This result is a consequence of a general result for Malcev-Neumann series [18, Theorem 
3.1.7]. As a consequence, the series expansion of 1/(1 — /) for ord(/) > ord(l) is just f -\- 

-I- • • • . More generally, for two iterated Laurent series A and B with ord(A) < ord(i?), the 
expansion oll/{A — B) is 



In the field K{{xi, . . .Xn)), we define a total ordering on the variables, which produces a 
total ordering on its group of monomials. This total ordering plays a central role in series 
expansion. By thinking of iterated Laurent series as numbers, ord(/) > ord(l) means that / 
is much smaller than 1, or / = o(l). Similarly ord{B) > ord(^) means that B is much smaller 
than A, OT B ^ o{A). 



oo 



i=0 




For instance, in K{{xi,X2,Xs)), we have 
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The analogous situation for complex variables would be informally written as 1 y> Xi > 
> ■ ■ ■ ^ Xn when expanding rational functions into Laurent series, where :» means ''much 
greater". See [16] and [14, p. 231]. 

The the following three computational rules are frequently used in constant term evaluations. 
Let F,G e K{{xi, . . . ,Xn)). 

1. Linearity: CT-r.(aF + bG) — aCT-j. F + feCTj,. G, if a and b are independent of Xi. 

2. If F can be written as X]fe>o ^k^i, then CT F — F\^.^q. 

Xi 

d 

3. Res —F = 0. 

Xi OXi 

Remark 2.4. Depending on the working field, rational functions Q{xi,X2, ■ ■ ■ , Xm) may have 
as many as m! different expansions. More precisely, if cr is a permutation of [m], then (5(x) 
will have a unique expansion in K{{xai,Xa2, ■ ■ ■ , x^Jj)- The expansions of Q{'x) for different cr 
are usually different. So we need to specify the working field whenever a reciprocal comes into 
account. 

Iterated Laurent series is to obtained by defining a total ordering on its variables (this idea 
is not new, e.g., [14; 16]). In fact, it is a special kind of Malcev-Neumann series, which has 
been studied in [18], and has apphcations to MacMahon's partition analysis. 



3. MacMahon's Partition Analysis 

MacMahon's Partition Analysis is used for counting the solutions to a system of hnear 
Diophantine equations and inequalities, and the number of lattice points in a convex polytope. 
Such problems can be converted into evaluating the constant terms of certain Elliott-rational 
functions. This conversion has been known as MacMahon's partition analysis, and has been 
given a new life by Andrews et al. in a series of papers [1-9] . 

Definition 3.1. An Elliott-rational function is a rational function that can be written in such 
a way that its denominator can be factored into the products of one monomial minus another, 
with the monomial allowed. 

In the one- variable case, this concept reduces to the generating function of a quasi-polynomial. 

MacMahon's idea was to introduce new variables Ai, A2, • • . to replace linear constraints. For 
example, suppose we want to count the nonnegative integral solutions to the linear equation 
2ai — 3a2 -|- 03 -|- 2 = 0. We can compute the generating function of such solutions as the 
following: 
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Eai a2 as _ \ ^ prp \2ai-3a2+a3+2 ai ag^as 
Xi X2 ^3 — / ^ '-^ J- ^ -i-l -^3 ■ 



<»i.a2,a3>0 ai,a2,a3>0 

2ai-3a2+a3+2=0 

Now apply the formula for the sum of a geometric series. It becomes 

CT ^ 



A (1-A2xi)(l-A-%2)(1-Ax3)' 

The above expression is a power series in Xi but not in A. 

It is clear that if there are r linear equations, we can compute their solutions by introducing 
r variables A, short for Ai, . . . , A^. Thus counting sohitions of a system of linear Diophantine 
equations can be converted into evaluating the constant term of an Elliott-rational function. 

Theorem 3.2. If F is Elliott-rational, then the constant terms of F are still Elliott-rational. 

This result follows from "The method of Elliott" (see [13, p. 111-114]) developed from the 
following identity. Note that we have not specified the working field yet. 

Lemma 3.3 (Elliott Reduction Identity). For positive integers j andk, 

1 1/1 1 



(1 - xX^){l - yX-'') 1 - xyX^-^ \1 - xX^ 1 - yX'^ 
Elliott's argument is that after finitely many apphcations of the above identity to an Elliott- 
rational function, we will get a sum of rational functions, in which every denominator has 
either all factors of the form 1 — xA\ or all factors of the form 1 — y/A*. Now taking the 
constant term of each summand is easy. 

Theorem 3.2 reduces the evaluation of CTa F to the univariate case CT;^ F by iteration. 
Unfortunately, the Elliott reduction algorithm is not efficient in practice. Other algorithms 
have been developed, and computer programs have been set up, such as the "Omega" package 
[6] . But we can do much better by the partial fraction method and working in a field of iterated 
Laurent series. 

Before going further, let us review some of the work in [6]. The key ingredient in their 
argument is MacMahon's Omega operator 0>, which is defined by: 



oo 



5 E ■■■ E ^i,...,3.Ai^---A^^:=E---E^-.-.-' 

Sl= — oo Sr = — oo Sl=0 Sr=0 

where the domain of the Asi,...,sr is the field of rational functions over C in several complex 
variables and Aj are restricted to a neighborhood of the circle |Aj| = 1. In addition, the >lsi,...,s^ 
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are required to be such that any of the T" — 1 sums 

oo oo 
Sij =0 Si. =0 

is absolute convergent within the domain of the definition of ^si,...,sr- 
Another operator Q= is given by 



Sl = — OO 5r = — oo 



Andrews et al. emphasized in [6] that it is essential to treat everything analytically rather 
than formally because the method relies on unique Laurent series representations of rational 
functions. 

It is not hard to see that their definition always works if we are working in a ring such as the 
ring of formal power series in x with coefficients Laurent polynomials in A, where x is short 

for and A is short for Ai, . . . , A^. In fact, this approach was used by Han in [12]. 

By Theorem 3.2, it suffices to consider the case of r = 1, since the general case can be done 
by iteration. In the previous work by Andrews et al. and by Han, the problem was reduced 
to evaluating the constant term (with respect to A) of a rational function of the form 

^^■^^ ni<.<ji-A^-'x,)ni<,<ji-?/./A'=o' 

This treatment has assumed the obvious geometric expansion: 

oo oo 

--— = J]A^X and 

In other words, for each factor / in the denominator, / has positive powers in A indicates 
that the series expansion of 1// contains only nonnegative powers in A; and / has negative 
powers in A indicates that 1// contains only nonpositive powers in A. In our approach, these 
indications are dropped off after defining a total ordering. 

We find it better to do this kind of work in a certain field of iterated Laurent series, because 
in such a field, we can use the theory of partial fraction decompositions in K{X) for any field 
K and any variable A. We illustrate this idea by solving a problem in [6, p. 252] with the 
partial fraction method. 

Problem. Find all nonnegative integer solutions a, 6 to the inequahty 2a > 36. 
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Solution. First of all, using geometric series summations wc translate the problem into a form 
which MacMahon calls the crude generating function, namely 

1 



fix, y) := -^Y. ^'""'"^"y' = 



a,6>0,2a-36>0 " o,fe>0 " ^ 

where everything is regarded as a power series in x and y but not in A. 
Now by converting into partial fractions in A, we have 

1 _?/(! + Ax^y + >?x) ^ 1 + \x^y 



(1- A2x)(l- A-3y) (i-x3y2)(A3_y) 

Where the right-hand side of the above equation is expanded as a power series in x and y, the 

second term contains only nonnegative powers in A, and the first term, 

y{\ + \x^y + }?x) _ y A~^ + \~'^x^y + X'^x 
(1 — x^y'^){a^ — y) 1 — x^y"^ 1 — X~^y 

contains only negative powers in A. Thus by setting A = 1 in the second term, we obtain 

jy^'^y) (i-xv)(i-^)' 

By a geometric series expansion, it is easy to deduce that 

{ (a, 6) e : 2a > 36 } = { (m + n + \n/2\ , n) : (m, n) e }. 

□ 

In solving the above problem, we see that partial fraction decomposition helps in evaluating 
constant terms, and that only part of the partial fraction is needed. 

4. Algorithm by Partial Fraction Decomposition 

Working in the field of iterated Laurent series has two advantages. First, the expansion of a 
rational function into Laurent series is determined by the total ordering " on its monomials, 
so we can temporarily forget its expansion as long as we work in this field. Second, the fact 
that F is a rational function in A with coefficients in a certain field permits us to apply the 
theory of partial fraction decompositions. 

Note that the idea of using partial fraction decompositions in this context was first adopted 
by Stanley in [14, p. 229-231], but without the use of computers, this idea was thought to be 
impractical. 

MacMahon's partition analysis always works in a ring like K[K, A~^][[x]], where is short 
for A^^, . . . , A~^. This ring can be embedded into a field of iterated Laurent series, such as 
i^((A,x)). 



A FAST ALGORITHM FOR MACMAHON'S PARTITION ANALYSIS 9 

While working in the field of iterated Laurent series, it is convenient to use the operator 
PT;^, which is formally defined by 

oo oo 

FT J2 anA" = J]a„A", 

n=— oo n=0 

whose validity is justified by the fundamental structure theorem. 
MacMahon's operators can be realized as the following. 



(4.1) OF(A,x) = PTF(A,x) 



A=(l,.-,1) 



(4.2) l]F(A,x) =CTF(A,x) = PTF(A,x 

A A A=(0,...,0) 

So it suffices to find PTa-F. In fact, it is well-known that Q> can be realized by Q= by 
introducing new variables, just as the PT operators can be reahzed by the CT operators (see 
[18, Ch. 1]). So either an algorithm for PTa-F or an algorithm for CTa-F will be sufficient 
for our purpose. Generally speaking, PT is more suitable for the algorithm, and CT is more 
suitable for theoretical analysis. 

Now we need an algorithm to evaluate PT^ F{\) with 

F(X) = ^(^) 

^ ^ ni<.<„(A^-^ - z.) 

where -P(A) is a polynomial in A, ji are nonnegative integers, and Zi are independent of A. Note 
that we allow Zi to be zero, so that the case of P{X) being a Laurent polynomial is covered. 
Our approach is different from the previous algorithms, which deal with rational functions 
expressed as in (3.1) (the difference will be further explained in the next section). It based 
on the following known fact, which says that once the partial fraction decomposition of F is 
given, PTx F can be read off immediately. 

Theorem 4.1. Suppose that the factors in the denominator of F are pairwise relatively prime, 
and that the partial fraction decomposition of F is 



where /(A) is a polynomial in X, and Pi{X) is a polynomial of degree less than ji for each i. 
Then 

(4.3) PTF = /(A) + E#^' 

i 

where the sum ranges over all i such that Zi -< X^^ . 



10 GUOCE XIN 

Proof. The condition that Zi is independent of A imphes that either A-'* -< Zi or Zi -< V\ In 
the former case, we observe that the expansion of pi{X)/{V^ — Zi) into Laurent series contains 
only negative powers in A, hence has no contribution when applying PTa- In the latter case, 
the expansion contains only nonnegative powers in A. Thus the the theorem follows. □ 

To apply Theorem 4.1, we need to know the partial fraction decompositions of the given 
rational function. In fact, we need only part of the partial fraction decompositions. Thus we 
need an efficient algorithm for the partial fraction decompositions. More ideally, an algorithm 
that only give us the necessary parts. The classical algorithm does not seem to work nicely. 
We use the new algorithm in [17] developed from the following Theorem 4.2. 

To state the theorem, we need some concepts. Let i^T be a field. For N,D & K[t] with D ^ 0, 
N/ D can be uniquely written as the summation of a polynomial p and a proper fraction (or 
rational function) r/D. We denote by Poly(A^/D) the polynomial part, which is p, and by 
Frac{N/D) the fractional part, which is r/D. 

Suppose that N,D e K[t\ and D is factored into pairwise relatively prime factors D — 
Di - ■ ■ Dk. Then the ppfraction (short for polynomial and proper fraction) expansion of N/D 
with respect to Di, . . . , is the decomposition of N/D as 

N/D^p + ri/D, + ---rk/Dk 

such that p,ri are polynomials and deg(rj) < deg(Dj) for every i. We denote the above ri/D^ 
by Frac{N/D, Di), the fractional part of N/D with respect to D^. 

Theorem 4.2 (Theorem 2.3 [17]). For any N,D e K[t] with D ^ 0, if Di, . . . Dk e K[t] are 
pairwise relatively prime, and D = Di - ■ ■ Dk, then 

N ^ , fN\ ^ fN ^\ ^ fN ^\ 

- = Poly (- j + F«c (-. A j + ■ ■ ■ + f«c (-, j 

is the ppfraction expansion of N/D with respect to {Di, . . . , D^). Moreover, ifl/{DiDj) — 
Si/Di +pi/Di, then 

Frac{N/D, D^) = Frac{Ns2S^ ■ ■ ■ Sk/D^). 

By Theorem 4.2, we need two formulas to develop our algorithm. One is for the fractional 
part of p(A) / (A-' — a), and the other for the partial fraction decomposition of (A-' — a)'~^(A'^— 6)~^. 
These are given as Propositions 4.3 and 4.6 respectively. 

Let n mod k be the remainder of n when divided by k. We have 
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Proposition 4.3. The fracAAonal part of p{\)/{V — a) can he obtained by replacing X'^ with 
^{dvaod j)^\d/j\ ^^^-j dividing the result by — a. 

Proof. By linearity, it suffice to sfiow tfiat tlie remainder of A'^ wfien divided by A-' — a equals 
^{d mod j)^[d/j\ ^ ^jjjpjj jg trivial. □ 

It is easy to sec that this operation takes time linear in the number of nonzero terms of p(A), 
where we assume fast arithmetic operations. 

Remark 4.4. Observe that the numerator of the fractional part of p(A) / (A-' — a) is always a 
Laurent polynomial in all variables. 

Lemma 4.5. For positive integers j and k, if a'' ^ V , then the following is a partial fraction 
expansion. 

^ ' {Xi-a){X^-h) y-a^ \ X^-h ) W-a^ \ X^ - a 

Proof. First we show that if ^ , then A-' — a and X^ — b are relatively prime. If not, 
say ^ is their common root in a field extension, then = a and = b. Thus we have 
^ (^^j-^k ^ ^jk ^ ^^ky ^]jp contradiction. 

We have 

h> -a^ _ X^^ - a'' X^'' - h> 

{Xi - a){X>' - b) ~ {X^ -a){X^-h) ~ {XJ - a){X'' - b) 

~ X^-b Xi -a ■ 

Now the polynomial part of (^j^nXAfe-b) clearly 0. Thus the sum of the polynomial parts of 
the two terms on the right side of the above equation also equals 0. So taking the fractional 
part of both sides and then dividing both sides by V — gives the desired result. □ 

Now if gcd(j, k) is not 1, then we can replace A^'^'^^^'*^-' with and apply the above lemma. 
This gives us the following result. 
Let 

J^{X^-a,X'-b)^ ^^=\, r, , 

where / = j/ gcd(j, k) and k' — k/ gcd(j, k). 

Proposition 4.6. For positive integers j andk, if a^ ^ W , then we have 

(4.5) Frsc f,, ' X'-a)= F.c (^(Vp^) , 

^ ' \{Xi - a){X^ - b) J \ X^-a J 
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Remark 4.7. Note that a similar result appeared in [6, Theorem 1], but their proof was 
lengthy. 

Now by Theorem 4.2, we have the following: 

Theorem 4.8. With the notation of Theorem 4-1, the polynomial Ps{X) equals the remainder 
of 

n 

P(A) n H^''-as,X''-ai), 

when divided by X^* — zi as a polynomial in X. 

In Theorem 4.1, we assumed that A''* — Zi and A-'* — z^ are relatively prime. Now let us 
consider the case that A''* — Zi and A-'* — Zk have a nontrivial common factor. This happens 
if and only if zf" = zl% which can be easily checked. Andrews et al. [6] suggested that we 
temporarily regard Zi and zj as two different variables. After the computation, we replace 
them. We find an alternate approach, which has been implemented in our computer program 
and will be discussed in the next section. 

Thus the above argument. Theorem 4.1, and Theorem 4.8 together will give us an fast 
algorithm for evaluating CT^ F. 

Remark 4.9. Prom Remark 4.4, Theorem 4.1, and Theorem 4.8, we see that PT^ F is EUiott- 
rational when F is. This is another way to prove Theorem 3.2. 

Example 4.10. Count all triples (a, 6, c) in such that they satisfy the triangle inequalities. 

Similar problems have been done, such as counting non-congruent triangles with integral 
side lengths [15, Exercise 4.16]. We are going to illustrate our new approach by this example. 

Solution. We solve the following three Diophantine inequalities: a + b — c>0, b + c — a>0, 
and c + a — b > 0. The generating function of these solutions is equal to fl^F{A, x), where 

^^^^ " ((1 - AiA3Xi/A2)(l - AiA2X2/A3)(l - A2A3X3/A1)) ■ 

Although F{A, x) is in K[A, A~^] [[x]], we shall work in a field of iterated Laurent series. We 
will chose K{{A,x.)). and K{{A,Xs,X2,Xi)) and compare the results. 

We will apply fl> to A3, A2, Ai subsequently. The first step, applying Q> to A3 makes no 
difference for the two working fields. Applying Theorems 4.1 and 4.8 to the factors of F{A) 
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containing A3, wc get 
fi>,A3i^(A,x) = 

{Xi^xi - Xlxs) (1 - Xi^X2Xi) {Xixi - A2) (Ai^a;i - X2^X3) (l - X2^X2X3) (Ai - X2X3) 

Denote by Fi and F2 the above two summands. At this stage, we note that the expansion 
of [Xi^xi — Al^a) ^ does not exist in ir[A, A~^][[x]], and generally there is no advantage in 
getting rid of the factor Afxi — A|x3 in the denominator by combining the above two summands 
into one rational function. This will be further explained in the next section. 

Now when applying il> on A2 to Fi and F2, the results are different for the two working fields. 
Let us look at Fi, especially the expansion of (Ai^xi — A2a;3) ^. The expansion in K{{A,x.)) 
contains only nonnegative powers in A2, while the expansion in K{{A,X3,X2,Xi)) contains only 
negative powers in A2. The situation for F2 is similar. The conclusion is that working in 
K{{A,X3,X2,Xi)) is better: applying f]> on A2 to Fi will gives us 0, and we have 



Q>,A3,A2i^(A,x) 



Ai (Ai + 0:3) X2 _^ X1X3 



{-X3 + Xi^X2) (-1 + Ai^a;2a;i) {X2X3 - 1) [-X3 + Xi^X2) (-1 + X1X3) {-X3 + Ai) 
Applying Q> on Ai to the two summands of the above equation and simplifying gives us 

X3X2XI 

□ 

Remark 4.11. It is left to the reader to check that for the above example, C((A, X2, X3, xi)) 
is the best working field. 

5. The Maple Package 

Lemma 5.1. Let ji be positive integers and let Zi he monomials. If X^^ — Z\ is not relatively 
prime to X^^ — Z2, nor to X^^ — Z3, then X^"^ — Z2 and A^* — Z3 are not relatively prime. 

Proof. By the proof of Lemma 4.5, we have z''^ = and z^^ = z^^ . It is easy to see that 
Z2^-^^ = z^^-'^ . Now the fact that Zi is a monomial (and hence has coefficient 1) implies that 

4' = 4'- □ 



14 GUOCE XIN 

Thus to obtain the complete algorithm, we need to handle the situation when A-'^ — Zi, . . . , 
X^*' — Zk are not relatively prime to each other. For this situation, we have not succeeded in 
applying the suggestion of the last section: we tried to let Zi = ZiVi and do the computation, 
and finally replace Vi with 1. But the problem is that the last step can only be done after 
simplification, for which the rational function will be too big for Maple to deal with. The 
following example explains why this is not a fast approach: evaluating 

1 



(l-Ax)io(l-y/A)8' 

Our current program uses a modified ppfraction expansion as follows. Suppose that N, D, 
Pi belong to K[X], and that P — pi- • - pk is relatively prime to D. Then we can obtain a 
formula for Frac(iV/PD, P) satisfying our needs: 

Write N/{piD) = ri/pi + Ni/D with deg(ri) < deg(pi). Then ri/pi = frac{N/piD,pi) can 
be easily obtained. 

Now write A^i/(p2^) = + N2/D. Then N/{piP2D) = ri/(piP2) + + N2/D. 
In general, we have 

I ' ' ' I ~r , 



Pi---PkD pi---pk Pk D 

with deg(ri) < deg(pi). Now it is easy to see that 

\PD J Pi---Pk Pk 
The recurrence formula for and Ni is given by 

n ^ fN,^, \ Ni^i-nD 

— — Frac — 7r,Pi and Ni 



Pi \PiD ) Pi 

where A^o = ■ Note that we shall let Maple compute Ni with respect to A. 
Now we can give the algorithm for computing PT^ F{X) as follows. 

(1) Collect the factors in the denominator of F into several groups, such that the factors 
in different groups are relatively prime and factors in a same group are not. 

(2) For each group having a contribution, find its corresponding fractional part of F . 

(3) Take the sum of the results obtained from step 2, and add the polynomial part of F . 

Remark 5.2. We will simplify only if needed. 

The factors in the denominator of F that are independent of A should be factored out to 
speed up the calculation. This has been implemented in our computer program. 

The algorithm for OaF(A, x) is described as follows. 
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(1) Fix a total ordering on x and a total ordering on A. Suppose we are working in C((A, x)). 

(2) Eliminate Ai by computing PT^^ F and then replacing Ai with 1. 

(3) For each rational functions obtained from step 2, eliminate A2. 

(4) Eliminate all the A's, and finally simplify. 

This approach partially solves the "run-time explosion" problem existing in Omega Calculus. 
Let us analyze a simple situation by considering Q>F(A), where 

^ ^ nr=i(i-VA)n;=i(i-%A)- 

The result after eliminating A and combining terms will have a denominator of mn factors: 
(1 — XiUj) with 1 < i < m and 1 < j < n. Such factors potentially contain the other variables 
that are going to be eliminated. This explains the existence of the run-time explosion problem. 

In our approach, the result after eliminating A will be a sum of n rational functions (with 
a possible polynomial part), each with a denominator of m + n factors. Now it is crucial that 
for each rational function, we can apply the theory of iterated Laurent series to eliminate the 
other variables. 

A Maple package implementing the above algorithm is available onhne at [19]. Here is an 
example of how to use this program after downloading this package. The current program uses 
E_Oge(F, X, A) to compute Q> F(A, x) in the field C((A, x)), where x is reahzed by [xi, • • • , x„] 
in maple and A is realized similarly. 

Example 5.3. Compute the generating function of k-gon partitions, which are partitions 
that can be the side lengths of a k-gon. 

This problem was first studied in [4], where the generating functions of k-gon partitions are 
obtained only for /c < 6 by using the authors' Omega package. We will discuss in the next 
section about their formula for general k. 

In the following F{k) is the crude generating function of k-gon partitions 

F{k) "^'"^^ 



(1 - x,fj{l - x,^) •••(!- - x,^) ' 

where we use Oj to replace Aj. The function test(/c) computes Jl> F{k) and gives its normal 
expression. 

> read "Ell.mpl"; 

> F:=proc(k) 

> product(l-q*a[k]*a[i-l]/a[i] ,i=2. .k-1); 
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> q/a [1] / ( ( l-q*a [k] /a [1] ) *%* ( l-q*a [k- 1] /a [k] ) ) ; 

> end: 

> va:=proc(k) seq(a[i] , i=l . .k) end: 

> F(3); 

V «1 / V 02 / V «3 / 

> E_Oge(%, [q] , [va(3)]) ; 

> test : =proc (n) F (n) ; va(n) ; E_Oge (%% , [q] , [%] ) ; normal (%) ; end : 

> test (3); 



(^4 _ 1) (^3 _ 1) (^2 _ 1) 

> test (4); 

{q' -q' + l) 



+ (52_l)2(^3 + ^2 + ^+l)(^2_^+l) 

> test (5); 

+ + + q'^ + + ^-q'^ + q^ + q'^ + q + l)q^ 
~ (g3 - 1) (-1 + q)(q^-q3 + q- 1) (g6 + ^5 _ ^ _ i) (_i + ^8) + (^2 + 

> test (6); 



> time (test (7) ) ; 

34.765 

All of the above are done in a personal computer. The time of test (7) is measured by 
seconds. 

Example 5.4. A Putnam problem (B3 on the 2000 Putnam examination) was generalized in 
[7]. The generalized problems are converted to evaluating CP{T, k, c) = fl> P{T, k, c) where 

^' ^ (1 - xi(ai . • • arra;^'^^-'^^^^) •••(!- xr{a^ ■ ■ ■ arra^^'^^-'^^^^) 
for A; > c and for A; < c we have 

P{T, k, c) = — - 7, (fc(T-i)+c)x" 

[1 — x\[ai ■ ■ ■ ax) Oi • • • (1 — XT(ai • • -ot) ) 

Note that the case of A; = c is trivial. 
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> read "Ell.mpl"; 

> P:=proc(T,k,c) if k>c then 

> 1/product (1-x [i] *product (a [ j ] "k , j=l . . T) /a [i] ^ (k* (T-1) +c) , i=l . . T) 

> else 

> 1/product (1-x [i] /product (a [j]"k,j=l. .T)*a[i] "(k+CT-D+c) ,1=1- -T) 

> end: 

> va:=proc(T) seq(a[i] , 1=1 . .T) end: 

> vx:=proc(T) seq(x [i] , 1=1 . .T) end: 

> CP:=proc(T,k,c) E_Oge(P(T,k,c) , [vx(T)] , [va(T)] ); normal (%) ; end; 

The case T — 3, k — 2 and c = 1 is given explicitly. 



> P(3,2,l); 



> CP(3,2,1); 



Xl^X2^X3^ + Xi'x2'x3' + X2XsX\ + XiX\X2, + 1 



(-1 + X-iX<i^X^) (-1 + X2X-^X^) {x<i^X-i^X\ — 1) 

> time(CP(3,3,l)) ; 



> time(CP(3,2,3)) ; 



> time(CP(3,l,3)) ; 



> time(CP(3,l,4)) ; 



> time(CP(3,l,5)) ; 



0.171 



0.391 



0.235 



0.686 



2.172 



> CP(4, 1,3) : factor (7.) ; 
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Table 1 . Comparison of our new method and the Omega package 



run-time for CP 


(3,3,1) 


(3,2,3) 


(3,1,3) 


(3,1,4) 


(3,1,5) 


(4,1,3) 


New inotliod 


0.171 


0.391 


0.235 


0.686 


2.172 


11.110 


Omega package 


7.14 


58.95 


100.55 


643.86 







Maple will give us the following result: 

(a;4^X3^a;2^a;i^ + 1 + X2X4X3X1) 

{X4X3XI^X2 — 1) {XiX4X3^X2 — 1) {XiX4^X3X2 — 1) {x2^X4X3Xi — 1) 

X2 X4 X:i Xi + X2 X4 X^ Xi + X4 Xs X2 Xi + X2 X4 X^ Xl + X2 X4 X:i Xi 

+ X\X4X-iX2 + X\X4X'iX2 + X2X4X-iX\' + X4X-iX2X\ + X\X4X'iX2 

+ X2X4XzX\ + X4X-iX2X\ + X-^X^XzXx + X4X-iX2X\ + 3:2^X4X3X1 + l) 

It will take Maple more than 5 minutes to evaluate CP{4:, 2, 3). 

We give a detailed comparison of the new algorithm and the Omega package in Table 1, 
where the unit of the run-time is seconds. Note that programs are not running on the same 
computer, and that the data for the Omega package comes from [7]. 

Example 5.5. Magic squares of order n are n by n matrices with integral entries such that 
all the row sums and column sums are equal. 

The crude generating function for Magic square is: 

t TTTT t 

1 - t{Xi ■ ■ ■ XnUl ■■■iJin) H" 1 - Xij/{XiHj) 

t J- J 1 

When evaluating the constant term in A and /.t's, we use E_Oeq instead of E_Oge. Our maple 
program will reproduce the result for n = 3 quickly. For n = 4, we get a sum of 96 simple 
rational functions, which is less than the 256 of the Omega package. Moreover, if we set 
Xjj = X at the beginning, and finally replace x with 1, our program will reproduce the formula 
for the case n = 5 in about two minutes, which is the generating function of the row sums for 
magic squares of order 5 [15, p. 234]. 

6. Ways to Accelerate the Program 

Our program should have been accelerated by several ways, which are not implemented 
due to the author's lack of programming skills. These ways are list as follows and explained 
by examples, we will manage to reduce of the number of rational functions of the output. 
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since the simplification of a sum of many rational functions is a bottle neck for Maple (also 
Matliematica) . 

(1) The order of the variables to eliminate can make a difference for the computational 
time. 

(2) The total ordering on the x's can make a difference for the computational time, as will 
the total ordering on the A's. 

(3) The following alternative formula of (4.3) can simphfy the computation: 



where the sum ranges over all i such that Zi >- 

(1) is a well-known fact. To take advantage from it, we use the fact that the number of 
rational functions produced by eliminating Aj is equal to the number c/(Aj) of factors in the 
denominator of F that have contributions with respect to Aj. If cf{Xig) is the smallest among 
all the cf{\i), then we shall ehminate the Ai^ first. Note that This way does not guarantee the 
best result. 

The first part of (2) can be explained by Example 4.10, which gives a simple example of how 
to take advantage of (2). The exact description will take time. The second part is similar [18, 
Example 2.5.13]. 

Using (3) might produce fewer rational functions. This happens when the denominator of F 
has more factors with contribution than those factors without contribution. See the following 
example. 

Example 6.1. Count all /c-gons with nonnegative integral side lengths, which are not required 
to be in an increasing order. 

Solution. Suppose the side lengths of a /c-gon is given by ai, . . . , a^. Then we have k inequali- 
ties, Oi -|- • — h Ojfc > 2aj for all i. 

Using formula (6.1) we can compute the generating function of /c-gons without computer. 
The eliminating order is A^, . . . , Ai. 



(6.1) 





1 



(1 - xiAi • • • Ajt/Af ) • • • (1 - XfeAi • • • Afc/A^) 



1 



(1 - xiAi • • • Ajfc_i/Af) • • • (1 - Xk-i\i ■ • • Afe_i/A^_J(l - XfeAi • • • Afe_i) 
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XkXi ■ ■ ■ Afc-i 



(1 - xiXkXl ■ ■ ■ Xl_^/ Af ) • • • (1 - Xk-iXkXl ■•■Xl_^/ A2_ J (1 - XfeAi • • • Ajk_i) ' 



Now notice that Jl> acting on the second term is simply obtained by replacing Aj with 1. 
Repeat the above computation, we get the final generating function: 



However, it is probably better not to use (6.1) if the total degree of those factors without a 
contribution is much greater than those factors with a contribution. Also note that this formula 
is not easy to apply for the CT operator, we shall use CTa-F(A) = CTaF(1/A) instead. 

There are also ways that may speed up the computation, but are not easy to implement by 
the computer. The following example is simplified by using different parameters for a given 
problem. 

Example 6.2. Count /c-gon partitions (revisited). 

An exact formula for the generating function of A:-gon partitions was given in [4, Theorem 
1]. Here we give a simple proof by using different parameters and formula (6.1). 

Solution. The problem is to find all (ai, . . . , a^) e P*^ such that 1 < ai < a2 < • • • < a*;, and 
ai H h flfc-i > flfc. 

Let hi = ai — 1, h2 = 02 — oi,. . . , hk = Ok — Ofc-i- Then = 1 + 6i + ■ ■ • + 6i for all i, and it 

suffices to find all hi such that hi > 0, and — 3 + (/c — 2)6i + (A; — 3)^2 H > bk- Thus 

the generating function for these hi are given by 




i=l 



□ 



> (1 - xiX^-^) • • • (1 - Xfe_2A)(l - Xfe_i)(l - Xk/X) 



(1 - xi) • • • (1 - Xk) (1 - a;ix^2)(l - X2X^3) • • • (1 - Xk-i){l - Xk) ' 



Now it is easy to convert this formula to [4, Theorem 1]. 



□ 



Acknowledgment. The author is very grateful to his advisor Ira Gessel. 



a fast algorithm for macmahon's partition analysis 21 

References 

1. G. E. Andrews, MacMahon's partition analysis. I. The lecture hall partition theorem, Math- 
ematical essays in honor of Gian-Carlo Rota (Cambridge, MA, 1996), Progr. Math., vol. 
161, Birkhauser Boston, Boston, MA, 1998, pp. 1-22. 

2. , MacMahon's partition analysis. II. Fundamental theorems, Ann. Comb. 4 (2000), 

327-338, Conference on Combinatorics and Physics (Los Alamos, NM, 1998). 

3. G. E. Andrews and P. Paule, MacMahon's partition analysis. IV. Hypergeometric mul- 
tisums, Sem. Lothar. Combin. 42 (1999), Art. B42i, 24 pp. (electronic). The Andrews 
Festschrift (Maratca, 1998). 

4. G. E. Andrews, P. Panic, and A. Ricsc, MacMahon's partition analysis. IX. k-gon parti- 
tions. Bull. Austral. Math. Soc. 64 (2001), 321-329. 

5. , MacMahon's partition analysis: the Omega package, European J. Combin. 22 

(2001). 

6. , MacMahon's partition analysis. VI. A new reduction algorithm, Ann. Comb. 5 

(2001), 251-270, Dedicated to the memory of Gian-Carlo Rota (Tianjin, 1999). 

7. , MacMahon's partition analysis. VII. Constrained compositions, ^-series with ap- 
plications to combinatorics, number theory, and physics (Urbana, IL, 2000), Contemp. 
Math., vol. 291, Amer. Math. Soc, Providence, RI, 2001, pp. 11-27. 

8. , MacMahon's partition analysis. VIII. Plane partition diamonds. Adv. in Appl. 

Math. 27 (2001), 231-242, Special issue in honor of Dominique Foata's 65th birthday 
(Philadelphia, PA, 2000). 

9. G. E. Andrews, P. Paulc, A. Riese, and V. Strchl, MacMahon's partition analysis. 
V. Bijections, recursions, and magic squares. Algebraic combinatorics and applications 
(GoBweinstein, 1999), Springer, Berlin, 2001, pp. 1-39. 

10. W. Baldoni-Silva and M. Vergne, Residues formulae for volumes and Ehrhart polynomials 
of convex polytopes, preprint, arXivrmath. CO/0103097, 2001. 

11. C. S. Chan, D. P. Robbins, and D. S. Yuen, On the volume of a certain polytope. Experi- 
ment. Math. 9 (2000), no. 1, 91-99. 

12. G. N. Han, A general algorithm for the MacMahon Omega operator, Ann. of Comb. 7 
(2003), 467-480. 

13. P. A. MacMahon, Combinatory Analysis, vol. 2, Cambridge University Press, Cambridge, 
1915-1916, Reprinted: Chelsea, New York, 1960. 

14. R. P. Stanley, Combinatorial reciprocity theorems. Adv. in Math. 14 (1974), 194-253. 



22 GUOCE XIN 

15. , Enumerative Combinatorics, 2 ed., vol. 1, Cambridge University Press, 1997. 

16. K. G. Wilson, Proof of a conjecture by Dyson, J. Math. Phys. 3 (1962), 1040-1043. 

17. G. Xin, A fast algorithm for partial fraction decompositions, J. Symbolic Comput. (2004), 
submitted. 

18. G. Xin, The Ring of Malcev- Neumann Series and The Residue Theorem, Ph.D. thesis, 
Brandeis University, 2004, arXiv:math.CO/0405133. 

19. Guoce Xin, A Maple package for the constant terms of Elliott-rational functions, 2004, 
http : / / www.br andeis .edu / ~ maxima / maple / Ell. mpl , . 

20. D. Zeilberger, Proof of a conjecture of Chan, Robbins, and Yuen, Elec. Trans. Numer. 
Anal. 9 (1999), 147-148. 



Department of Mathematics, Brandeis University, Waltham MA 02454-9110 
E-mail address: maxima@brandeis.edu 



